Anatomical harmonics basis based brain source localization with application to epilepsy

Brain Source Localization (BSL) using Electroencephalogram (EEG) has been a useful noninvasive modality for the diagnosis of epileptogenic zones, study of evoked related potentials, and brain disorders. The inverse solution of BSL is limited by high computational cost and localization error. The performance is additionally limited by head shape assumption and the corresponding harmonics basis function. In this work, an anatomical harmonics basis (Spherical Harmonics (SH), and more particularly Head Harmonics (H2)) based BSL is presented. The spatio-temporal four shell head model is formulated in SH and H2 domain. The anatomical harmonics domain formulation leads to dimensionality reduction and increased contribution of source eigenvalues, resulting in decreased computation and increased accuracy respectively. The performance of spatial subspace based Multiple Signal Classification (MUSIC) and Recursively Applied and Projected (RAP)-MUSIC method is compared with the proposed SH and H2 counterparts on simulated data. SH and H2 domain processing effectively resolves the problem of high computational cost without sacrificing the inverse source localization accuracy. The proposed H2 MUSIC was additionally validated for epileptogenic zone localization on clinical EEG data. The proposed framework offers an effective solution to clinicians in automated and time efficient seizure localization.

The use of an accurate and efficient algorithm for Brain Source Localization (BSL) using Electroencephalography (EEG) measurements has been used in various neuroscience applications such as pre-surgical mapping in patients undergoing resection of epileptogenic zones 1-4 , Brain Computer Interface (BCI) based prosthetic limbs 5,6 , and attention deficit hyperactive disorder 7,8 . Among these diseases, epilepsy is the most important and common neurological disorder as 1% of the world population is suffering from it. To solve this problem, EEG is regarded as the most commonly used non-invasive diagnosis tool due to its high temporal resolution, portability and cost-effectiveness. The process of brain source localization is carried out in two phases, which are forward modeling and inverse modeling.
Forward modeling involves estimation of the scalp potentials, given head model and current source. The volume conductor head model is chosen based on application, accuracy and computational complexity. The two kinds of volume conductor head models include numerical and analytical model. The numerical models are based on realistic modeling of human head with more complexity and low localization error. Boundary Element Method (BEM) 9 , Finite Element Method (FEM) 10 , and Finite Difference Method (FDM) 11 belong to numerical model category. In analytical modeling, head is approximated as a set of concentric spheres, each with homogeneous conductivity. This approximation results in low computational complexity. In literature, a single-shell head model 12 with a homogeneous conductivity was reported first. Later, three-shell 13 volume conductor head model was introduced where conductivity of skull was found to be lower than that of scalp/brain tissues. A more accurate four shell head model 14 with additional CerebroSpinal Fluid (CSF) layer between brain and skull, was introduced thereafter that forms the basis of the current work.
The recorded scalp potentials are further utilized to solve the inverse problem for active source localization. BSL methods are widely classified into distributed source (dipole-imaging) model and the dipole-fitting model. The distributed source model assumes that there is a large number of sources confined in an active region of the brain. The BSL, in this case, refers to the estimation of active source amplitudes and orientations using linear optimization techniques such as Bayesian methods, Minimum Norm Estimates (MNE), Weighted Minimum Norm Estimates (WMNE), Low Resolution Electrical Tomography (LORETA), and Local Auto Regressive

Four shell based EEG forward model in spatio-temporal domain
Spatio-temporal domain data model accounts for change in EEG potential w.r.t. space and time. Forward model is combination of head model and source model. The head model utilized herein is four layer concentric shell head model as shown in Fig. 1a, with the coordinate system in Fig. 1b. Brain, CSF, skull and scalp constitute the four concentric spherical layers with conductivity as σ 1 , σ 2 , σ 3 and σ 4 respectively. The brain comprises of homogeneous neural tissues and forms the innermost concentric sphere of radius br. This is surrounded by a concentric outer spherical shell of radius cr representing the CSF. This is followed by a third layer with radius dr for the skull. The outermost layer is the scalp of radius R. The BSL problem under consideration, utilizes Equivalent Current Dipole (ECD) source model where each dipole is parameterized by its location r p and dipole moment m p . I Sensors are placed over the scalp to record brain electrical activity due to P active dipole sources. The location of the ith sensor and pth source is given by respectively, where θ is elevation angle measured downward from positive Z axis, and φ is azimuth angle measured anticlockwise from positive X axis. All the dipole sources are assumed to be inside the brain compartment of the four shell model. The EEG potential V (r i , r p , t) measured at ith sensor location generated by pth active source dipole at a time instant t is given by 14 (1)  www.nature.com/scientificreports/ and (·) represents the vector dot product. Here {b, c, d} < 1 denote the relative radii of brain, CSF and skull with respect to the scalp radius. Eccentricity of the dipole is r p /R and f n (r p ) = (r p /R) n−1 . Radial unit vector is given by r o = r p /|r p | and t o represents tangential unit vector defined in terms of vector cross product as P n (cos �) is the legendre polynomial and P 1 n (cos �) is the associated legendre polynomial of order n ∈ {1, . . . , ∞} . The terms corresponding to n = 0 have zero contribution in (2) and is therefore not considered.
is the angle between ith sensor and pth dipole with cos given as Under the assumption of fixed dipole source location and orientation, dipole moment vector m p (t) can be written as m p (t) = e p s p (t) 16 , where e p = [e px e py e pz ] T is unit orientation vector and s p (t) is magnitude of the dipole moment. Substituting m p (t) in (2) and rearranging, the EEG potential at ith sensor due to pth dipole source can now be written as From (5), the total EEG potential at the ith electrode due to all the P active dipoles can be written as where g(r i , r p ) 1×3 is gain vector, expressed as for four shell model. For I electrodes and P dipoles, the EEG potential in (6), can be written in matrix form as Re-writing (8) in compact form, we have where elements of G(r p ) is given by (7), and [M] is referred as orientation moment matrix. In the presence of spatially and temporally zero mean white Gaussian noise Z with variance σ 2 , the four shell spatio-temporal data model in (9) for N s time snapshots can be written as where manifold matrix A = G M is collection of manifold vectors given by www.nature.com/scientificreports/ It is to be noted that matrix S is order independent and hence would remain same in the order limited scenario.

Manifold matrix
In this Section, the spatial domain manifold matrix is expressed in terms of spherical harmonics basis functions ( Y m n ) corresponding to sensor and source coordinates. Spherical harmonic addition theorem suggests 25 Substituting (14) in the gain vector g(r i , r p ) in (7), we have Radial and tangential orientation component of a dipole are expressed as respectively. For a radially oriented dipole, e pr = 1 and e pt = 0 . For a tangential oriented dipole e pt = 1 and e pr = 0 . A mixed oriented dipole has both radial and tangential component as non-zero. As the EEG reflects the electrical activity of the post-synaptic potentials of pyramidal neuron cells oriented perpendicular to the cortical surface, dipoles with radial orientation is considered herein. For a radially oriented dipole, the element a(r i , r p ) = g T (r i , r p )e p of a manifold vector can be written utilizing (15) and (16) as The expression in (17), can not be computed in practice due to the infinite series summation. EEG mode strength b n (r p ) determines the relative weighting of sensor and source spherical harmonics with the EEG signal order n. The parameter b n (r p ) is a function of head model and the source radial position and is plotted in Fig. 2 is defined by replacing the argument to p in (19). EEG mode strength B(r p ) is a diagonal matrix given by (12) a(r p ) = a(r 1 , r p ) a(r 2 , r p ) · · · a(r I , r p ) T = g T (r 1 , r p )e p g T (r 2 , r p )e p · · · g T (r I , r p )e p T = G(r p )e p (13) where α n = P 1 n (cos �) nP n (cos �) (15) e pr = r o · e p and e pt = t o · e p (16) www.nature.com/scientificreports/ The manifold matrix A I×P for a total of I electrodes and P dipoles can now be written as and , contain all sensors and dipoles position ( θ,φ)'s respectively. It may be noted from (21) that the manifold matrix is the product of SH basis functions corresponding to sensor coordinates, EEG mode strength and the SH basis functions corresponding to source coordinates.

Four shell based EEG forward model in spherical harmonics domain
In this Section, four shell based EEG forward model is presented in spherical harmonics domain. The spherical harmonic basis functions is first introduced followed by a detailed derivation of EEG signal decomposition in SH domain.
The spherical harmonics basis functions. The head model considered herein is four shell spherical model. Hence, the SH basis functions becomes a natural choice to a function defined on the head surface. In this work, SH basis functions are explored for forward problem formulation and subsequent localization, followed by a more accurate head harmonics 23,24 basis functions based approach. As EEG signals are real and discrete in nature, the real Spherical Fourier Transform (SFT) is applied instead of complex SFT. The real SFT pair for the discrete time EEG signal is given by 25 Here, dependency of V on (R, t) is omitted for notation simplicity. The sensor array order N a is dependent on spatial sampling scheme. Assuming nearly uniform sampling scheme with zero or negligible aliasing error, the number of sensors I should be at least κ(N a + 1) 2 , where κ assumes value in [1, 1.5] 26,27 . Equivalently, the order up to which the EEG potential can be recorded without aliasing is limited by number of sensors as It is to note that (26) presents an upper limit on array order that is function of number of sensors and the value of κ . The mathematical formulation of κ is beyond the scope of the present work.
The real valued spherical harmonics basis functions Y m n (θ, φ) of order n and degree m is defined as where, the normalization constant K m n is given by Spherical harmonics decomposition of EEG signal. In this Section, spatio-temporal forward model in (12) is reformulated in computationally efficient spherical harmonics domain. In practice, the continuous potential on the scalp surface is spatially sampled by placing the EEG electrodes as per the international electrode placement system. Hence, rewriting the forward SFT in (25) by replacing the integral, we have www.nature.com/scientificreports/ where the quadrature weights α nm i are chosen such that the approximation error involved in (29) is minimized. Substituting for V (θ i , φ i ) from inverse SFT (25) in (29), we have To ensure error-free sampling, following orthonormality condition must be satisfied.
Here δ denotes the Kronecker delta function. (31) can be rewritten in matrix form as and � a = (N a + 1) 2 . The quadrature weight matrix Q is estimated in the weighted least squares sense, given by A common choice for the weighting matrix W is an identity matrix. An approximated spherical fourier transform in (30) can now be written in matrix from as and the corresponding inverse SFT is given in matrix form as The quadrature weight matrix transforms the spatial domain forward data model to spherical harmonics domain data model. Multiplying Q from left to (12) results in SH domain data model as where the SH domain manifold matrix is written using (21) as A SH = QY(�)BY(�) . It is to be noted that the array manifold matrix of dimension I in spatio-temporal domain is reduced to dimension a , where a ≤ I . This reduction in dimensionality of the data is responsible for reduced computational cost in the transformed domain.

Four shell based EEG forward model in head harmonics domain
In this Section, four shell based EEG forward model is presented in head harmonics domain. The head harmonic basis functions is first introduced followed by a detailed derivation of EEG signal decomposition in H 2 domain.
The head harmonics basis function. The concept of spherical harmonics has been extended to develop a set of HemiSpherical Harmonics (HSH) basis functions that are defined over the unit hemisphere 28 . HSH basis functions have been found useful for representation of surface reflectance functions 29 , brain source localization 30 and surface description 31 . The HSH basis functions considers the EEG potential field to be defined on the hemisphere and therefore the human head is modeled as hemisphere 30 . However, the EEG signal acquired over head assumes the shape between hemisphere and sphere. Hence, based on realistic head dimension and EEG sensor placement system, another set of basis functions (called Head Harmonics (H 2 )) was proposed in 23,24 . It is to note that for EEG acquisition, the electrodes are placed over the head surface with an elevation angle in the range θ ∈ [0, 2π/3] range 27 . Therefore, it is more appropriate to utilize the H 2 basis functions (where θ ∈ [0, 2π/3] ) instead of SH basis functions (where θ ∈ [0, 2π] ). The real valued orthonormal H 2 basis functions ( H m n ) are expressed as 23 where K m n is a normalization constant and is given by www.nature.com/scientificreports/ and P m n (x) = P m n (1.33x − 0.33) is the shifted ALPs.
Head harmonics decomposition of EEG signal. The spatio-temporal forward model in (12), is reformulated in head harmonics domain. The SH decomposition of EEG signal in (38), assumes the potential field to be sampled on the sphere. However, the anatomical brain structure and the EEG sensor array placement suggest for the potential belonging to two third (i.e 120 • ) of the sphere only. Hence, the data model in (38)  It may be noted that the transformation matrix T H 2 = QH(�)Q transforms the spatial domain data model to the head harmonics counterpart. Additionally, the SH transformation matrix for converting the spatial domain data model to the spherical harmonics counterpart is T SH = QY H (�)Q = Q . In case of no transformation from the original spatial domain, the spatial transformation matrix T Spatial = I . The transformation matrices T SH and T H 2 , as detailed in Section 'Head Harmonics Decomposition of EEG Signal' , converts the spatio-temporal data model to SH and H 2 domain data model, respectively. In this process, the dimensionality of the data model is changed from I to � a = (N a + 1) 2 , where a ≤ I 26,27 . This reduction in dimensionality of the data is responsible for reduced computational cost in the transformed domain.

Spatial, SH and H 2 based inverse brain source localization algorithms
Following the development of four shell forward data model in χ ∈ Spatial, SH, H 2 domain, dipole fitting algorithms are presented herein for EEG inverse problem. In particular, subspace based MUSIC 16 and RAP-MUSIC 17 methods are formulated in the χ domain. The subspace based methods utilize the estimated noise or signal subspace to compute the corresponding cost functions. A peak in the cost function is attributed to the neural activity. It is to note that spatial MUSIC and RAP-MUSIC algorithms utilize the spatial domain data model in (12) while the proposed SH and H 2 based methods utilize the data model in (38) and (48) respectively. The spatial domain EEG signal V is first transformed to the desired χ domain counterpart by multiplying the corresponding transformation matrix T χ to obtain V χ = T χ V.
In the MUSIC algorithm, the estimates of active brain source location are obtained by projecting the χ domain array manifold vector a χ (r p ) = T χ a(r p ) over the estimated noise subspace. The χ domain MUSIC spectrum is expressed as www.nature.com/scientificreports/ where P χ is the noise subspace projection matrix obtained from eigen value decomposition of covariance matrix 16 . The MUSIC spectrum in (49) can also be written in a simplified manner as 23 where min (.) is the minimum eigenvalue of (.), and U χ is the Left Singular Vector (LSV) of T χ G(r p ) . One of the drawback of MUSIC algorithm is that it searches for multiple local peaks in the entire head volume, which makes the task time consuming and susceptible to false source localization. To overcome the limitation, RAP-MUSIC was proposed that extracts the global maxima at each recursion steps. For the kth recursion step, the χ domain RAP-MUSIC spectrum is defined as projects out the topography of the already found sources as with ⊥ A 0 = I . Matrix A k−1 concatenates the array manifold vectors of already found sources, and is given by

Simulated data analysis
Various experiments were conducted to illustrate the advantages of SH and H 2 domain based BSL algorithms. Both simulated and real EEG data were utilized for this purpose. In simulation, four layer concentric shell head model with radius of brain, CSF, skull, and scalp was considered to be 8.0 cm, 8.2 cm, 8.7 cm, and 9.2 cm respectively. The brain and scalp conductivities were set to 0.33 (�m) −1 . The conductivities of skull and CSFs were taken to be 1/80 and 5 times of that of the brain. A total of I = 128 Sensors were placed over head scalp utilizing 10 − 5% electrode placement system. The order N ref was set to 60 considering the diminishing nature of EEG mode strength for higher order. The number of discrete time samples N S was 200. The inter-grid gap was chosen to be 1 mm. The true source grid (i.e., the Region Of Interest (ROI)) was a 2-D axial slice at the depth of 3.2 cm from the top of the scalp surface as shown in Fig. 2b. The simulated sources were placed pseudo-randomly on ROI so that they were at least 3 cm from each other and at least 2 cm away from the center of ROI. Example of such a source distribution for P = 4 is shown in Fig. 2c. Active source location and orientation is assumed to be fixed. Potential data were generated considering the radial orientation of active brain sources. All experiments were conducted with L = 200 Monte-Carlo repetitions. It may be noted that the peak of the spectrum was explored only at ROI to reduce the computation load of the search over the entire volume. Spatial resolution. The capability of BSL algorithms in estimating the spatially closed active source location is analysed herein by plotting the respective cost function spectrum. Two dipole sources were placed at (2cm, −4cm, 6cm) and (4cm, −4cm, 6cm) on ROI. The inter source correlation and distance between the two sources was taken to be 0.2 and 2 cm respectively. The Signal to Noise Ratio (SNR) was set to 5dB. The order N a was taken to be 3 for both SH and H 2 domain. The 3-D view of spatial MUSIC, SH-MUSIC and H 2 -MUSIC cost function is plotted in Fig. 3 (a)-(c). It may be noted that the spatial MUSIC fails to localize closely spaced active dipoles at low SNR, providing a single false location. However, the proposed SH and H 2 counterparts localizes the two sources well with high resolution. This may be attributed to the increased contribution of the initial two eigenvalues of the covariance matrix corresponding to the two sources. The contribution of the initial two eigen  Additionally, the minimum distance between the two sources that have been localized via different BSL algorithms are quantified. Two sources were placed pseudo-randomly on ROI so that they were at least 2 cm away from origin. The inter source distance was varied from 0.5 cm to 2.0 cm. The minimum inter source distance is the shortest possible distance between two sources that can be localized. The minimum inter source distance and the corresponding localization error is illustrated in Fig. 4. It is to note that H 2 based MUSIC and RAP-MUSIC attains least localization error and minimum inter source distance when compared to spatial and SH counterparts. In particular, the H 2 RAP-MUSIC has the best performance. Localization error. In this Section, the Localization Error (LE) of spatial, SH and H 2 based BSL methods is presented. The LE for L iteration and P dipole source is defined as where (x, y) is the true source location and ( x, y) is the estimated source location. Simulation experiments were conducted to compute the LE with varying SNRs, array orders, number of active dipoles, and inter source correlation.
LE with varying SNR and array order for single source. The number of active dipole is taken to be P = 1 to study the effect of varying SNR and array order on LE. The source was placed in the pseudo-random manner. The LE is presented in Table 1 for MUSIC method in spatial, SH and H 2 domain. As expected, LE decreases consistently with increasing SNR. LE variation is additionally presented with varying array order N a ∈ {1, . . . , 10} in accordance with (26). High localization error of SH and H 2 MUSIC methods was observed for array order N a = 1 . This can be attributed to poor signal representation at lower order. The error further decreases with increase in order to give an optimum array order for a given number of sources. It is interesting to note that for single source, array  www.nature.com/scientificreports/ order N a = {2, 3} has the least localization error as highlighted in the Table 1 for SH/H 2 based localization. An increase in the localization error was observed thereafter thus giving an upper bound on the array order. This is consistent with the array signal processing 20 . In general, it may be observed that the H 2 MUSIC method outperforms the SH counterpart. For multiple source localization, the optimal array order with the least localization error is presented next.
LE with varying array order for two sources. Performance of spatial, SH and H 2 based MUSIC only method was explored in previous Section. In this Section, performance of MUSIC and RAP-MUSIC is presented for two active dipole sources with inter-source correlation as 0.2, placed pseudo randomly at 5 dB SNR. Localization error of MUSIC and RAP-MUSIC methods in the three domain with varying order N a is presented in Fig. 5a. As the spatial methods are independent of the array order, the localization error remains constant with order N a . It may be noted that H 2 , SH domain methods outperformed spatial counterparts. A valley in localization error plot may be observed at order 3 for the two active dipole sources. In all the three domains, the performance of RAP-MUSIC is better than MUSIC.

LE analysis for RAP-MUSIC with varying array order and number of active dipoles.
Optimum array order is experimentally estimated in this Section with varying number of active dipole sources. With an increase in the number of active sources, finding the local maxima in the MUSIC spectrum is susceptible to false source localization. Hence, LE results are presented in Table 2 for RAP MUSIC method only. The number of sources were varied from 1 to 5 at 5 dB SNR. The inter-source correlation between each adjacent pair was set to [0.20, 0.17, 0.17, 0.20]. It may be noted that the spatial domain RAP-MUSIC method localizes sources upto 3.
On the other hand, the proposed SH and H 2 methods successively localize all the five sources. An upward shift in optimal array order with an increase in the number of active sources was observed. This may be because of κ that appears to be dependent on the number of active dipole sources. Variation in the optimum array order with varying number of sources is highlighted in Table 2.  www.nature.com/scientificreports/ LE with varying inter source correlation. The subspace based source localization methods work on the prior assumption that the active sources are uncorrelated. When inter source correlation is significant (synchronous source case), BSL methods provide erroneous estimation of dipole sources. Therefore, to study the applicability of SH and H 2 based methods in the realistic scenario of synchronous source activation is of great significance. In the present simulation, the pearson correlation coefficient was varied from 0 to 0.8 for the case of two simultaneously active dipole sources at 5 dB SNR. The harmonics order for the SH and H 2 based methods was taken to be 3 as discussed in Fig. 5b. Effect of inter source correlation on the localization error is presented in Fig. 5b. It can be observed that the SH and H 2 domain based methods perform relatively better than their spatial domain counterparts in synchronous source localization. In particular the H 2 based method outperforms the rest.

Computation time. The computational efficiency of SH and H 2 domain based BSL methods such as MUSIC
and RAP-MUSIC is presented herein. The transformation matrices T SH and T H 2 , as detailed in Section 'Head Harmonics Decomposition of EEG Signal' , converts the spatio-temporal data model to SH and H 2 domain, respectively. In this process, the dimensionality of the data model is changed from I to � a = (N a + 1) 2 , where a ≤ I 26,27 . This reduction in dimensionality of the data is responsible for reduced computational cost in the transformed domain. The advantage of data model formulation in SH and H 2 domain in terms of dimensionality reduction is illustrated as a reduction in computation time. A 64-bit type system having Intel Core i7-6700 processor with 4 cores @ 3.40 GHz, and 16 GB RAM was utilized. Windows 10 Pro Version 1803 with MATLAB Version: R2019b was used for the evaluation. Computation time was measured using tic, toc MATLAB function at 5dB SNR. In the simulation, number of dipole sources were taken to be P = 2 . The average computation time was evaluated with varying sensor array order N a ∈ {1, . . . , 10} , as presented in Fig. 5c. It was observed that the computation time for SH and H 2 based methods is same. Present simulation demonstrate the effectiveness of SH and H 2 domain methods when compared with spatial domain methods. A significant reduction in the computation time for harmonics order N a ∈ [1, 5] is observed. It may be noted from Sect. 7.2 that the optimum harmonic order for two sources is 3. The high computational cost involved in solving the inverse problem using conventional spatial domain methods is a major obstacle in real-time EEG source localization.
Success rate. At a particular SNR value, the success rate of an algorithm is defined as the ratio of the number of localizations achieved with a 0 localization error to the total number of attempts. We further extend our performance evaluation by examining the success rate in % for spatial, SH and H 2 domain MUSIC methods. The number of Monte-Carlo trials were taken to be L = 1000. The order of SH and H 2 harmonics was set to 2. The success rate in % for a single source is presented in Fig. 6. At low SNR, significant difference in success rate can be observed between the SH, H 2 and spatial domain methods. In particular, the H 2 MUSIC is seen as the robust and most accurate source localization algorithm.

Real data analysis
For epilepsy surgery, localization of the seizure focus is time consuming and requires evaluation of neurological tests by experienced clinicians. In this Section, additional validation of the proposed head harmonic-based BSL is presented for epileptic seizure location using real clinical EEG data.
Subjects. An electrophysiological and imaging data of five patients with drug resistant epilepsy of long duration were included for study. The mean age of the cohort was 24.4 ± 7.33 with standard deviation range as 15 to 34 year. The patients were admitted to the Deenanath Mangeshkar Hospital & Research Center at Pune, Maharashtra, India. All of these patients were on multiple Anti-Seizure Medications (ASM). All patients underwent  EEG recordings. EEG recordings were done using standard 10-20 international system (19 electrodes along with ground and reference electrodes). Additional two electrodes T1 and T2 were added to detect anterior temporal discharges along with Electromyography (EMG) and Electrocardiography (ECG) electrodes. All the electrodes were pasted using standard gel and impedance was kept below 5k . Sampling frequency of EEG was set 2 KHz. Signals were passed through band pass filter of 1 Hz -70 Hz, followed by Notch filter to remove '50Hz' power line interference during recording. An epoch of 20 seconds was selected for source analysis by expert lab technician from the long duration recorded EEG, where first five seconds were pre-ictal recordings and another 15 seconds were seizure onset and propagation across the channels. Data pre-processing was performed offline on the selected epoch of data using the EEGLAB toolbox. The EMG and ECG channels present in the dataset were not considered.

Source localization.
Head harmonic based MUSIC algorithm is evaluated on the real epileptic seizure data as it provides all the peaks corresponding to active dipoles simultaneously. The data recorded from 21 channels is interpolated to 64 channel data using EEGLab toolbox 'spherical' interpolation method. The number of source was assumed to be two as eigen value decomposition of covariance matrix had 95% of its energy confined to a rank-2 subspace. Five epileptic patients data were analyzed. Epileptogenic zone was observed in the left temporal lobe for the patients P1, P2, P3 and P5. The patient P4 had epileptogenic zone in the right temporal lobe. Fig. 7 illustrates the H 2 MUSIC cost function plot of patient P1 having epiletogenic zone on the left temporal lobe. A pronounced cortical activation is focused in the corresponding epileptic region. The results are validated by the neurosurgeon who performed standard temporal lobe resection surgeries on left side in 4 patients (P1, P2, P3, P5) and on right side in 1 patient (P4). Our result demonstrates the effectiveness of H 2 domain algorithm in successfully localizing the true epileptic activation region on real EEG data.
Ethical approval. All methods were carried out in accordance with relevant guidelines and regulations approved by the institutional ethical committee of Deenanath Mangeshkar Hospital and Research Centre (Ref: DMHRC Code-IHR_2021_Apr_NK_403 , Dated: 2nd Apr 2021). A generic patient consent is taken but specific consent is waived off by ethics committee because this study did not alter the patient management and it's only involves studying of data acquired for other reason.

Conclusion
In this work, we proposed Spherical Harmonics (SH) and Head Harmonics (H 2 ) domain processing framework for active dipole localization. The four shell head model in spatio-temporal domain is formulated in computationally efficient SH and H 2 domain. A qualitative and quantitative performance comparison of Multiple Signal Classification (MUSIC) and Recursively Applied and Projected (RAP)-MUSIC methods in spatial, SH and H 2 domain is presented on simulated data. SH and H 2 domain processing effectively solves the problem of high computational cost without sacrificing the inverse source localization accuracy. The performance of H 2 based RAP-MUSIC was best when compared with SH and spatial domain counterparts. The performance dependency of SH and H 2 domain methods on array order was experimentally established. The proposed H 2 MUSIC was additionally validated for epileptogenic zone localization using clinical EEG data. The proposed framework offers an effective solution for clinicians in automated and time efficient seizure localization.